function [LMAP LE] = TurboDecode2(g, noisyDecIn,Lc, b_s, apriori, pos)
%b_s= block_size

x=noisyDecIn(1,pos);
%p1=noisyDecIn(2,:);
p2=noisyDecIn(3,:);

% generator polynomials parameters
[n,K] = size(g);
% number of memories
m = K-1;
% number of states
nb_st=2^m;

[trellis_matrix, parity_matrix]= trellis(g);
LMAP=zeros(b_s,1);
LE=zeros(b_s,1);


%Initialization Decoder MAP2
alpha= zeros(nb_st,b_s+1);
alpha(1,1)=1;
beta= zeros(nb_st,b_s+1);
beta(:,end)=1/(nb_st)*ones(nb_st,1);
    
gamma(1,1:b_s)= struct('gamma_matrix',zeros(size(trellis_matrix)));
gammaE(1,1:b_s)= struct('gamma_matrix',zeros(size(trellis_matrix)));

for k=1:b_s
    gamma(k).gamma_matrix= exp(1/2*(trellis_matrix.*(apriori(k) + Lc*x(k)') + Lc*parity_matrix.*p2(k)));
    gamma(k).gamma_matrix= gamma(k).gamma_matrix.*(trellis_matrix>0 | trellis_matrix<0);
%     gammaE(k).gamma_matrix= exp(1/2*Lc*parity_matrix*p2(k));
%     gammaE(k).gamma_matrix= gammaE(k).gamma_matrix.*(trellis_matrix>0 | trellis_matrix<0);
end

%Compute alpha1
for k=1:b_s
    alpha(:,k+1)= (gamma(k).gamma_matrix)*(alpha(:,k));
    alpha(:,k+1)=alpha(:,k+1)/sum(alpha(:,k+1));
end
%Normalize alpha1
% for i=1:length(alpha)
%     alpha(:,i)=alpha(:,i)/sum(alpha(:,i));
% end

%Compute beta1 with formula (15) 
for k=b_s:-1:1
    beta(:,k)= (gamma(k).gamma_matrix)'*(beta(:,k+1));
   beta(:,k)=beta(:,k)/sum(beta(:,k));
%    beta(:,k)=beta(:,k)/sum(alpha(:,k)'*(gamma(k).gamma_matrix)');
  
end
%Normalize beta1
% for i=1:length(beta)
%     beta(:,i)=beta(:,i)/sum(beta(:,i));
% end

for k=1:b_s
    sumSplus=alpha(:,k)'*(gamma(k).gamma_matrix.*(trellis_matrix>0))'*beta(:,k+1);
    sumSminus=alpha(:,k)'*(gamma(k).gamma_matrix.*(trellis_matrix<0))'*beta(:,k+1);
    LMAP(k)=log(sumSplus/sumSminus);
end

LE=LMAP-apriori-Lc*x';

end    
 